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Abstract 

Matheraatical models are an important tool for neuroscientists. During 
the last thirty years many papers have appeared on single neuron description 
and specifically on stochastic Integrate and Fire models. Analytical results 
have been proved and numerical and simulation methods have been developed 
for their study. Reviews appeared recently collect the main features of these 
models but do not focus on the methodologies employed to obtain them. Aim 
of this paper is to fill this gap by upgrading old reviews on this topic. The idea 
is to collect the existing methods and the available analytical results for the 
most common one dimensional stochastic Integrate and Fire models to make 
them available for studies on networks. An effort to unify the mathematical 
notations is also made. This review is divided in two parts: 

1. Derivation of the models with the list of the available closed forms 
expressions for their characterization; 

2. Presentation of the existing mathematical and statistical methods for 
the study of these models. 
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1 Introduction 

Progresses in experimental techniques, with the possibility to record simulta- 
neously from many neurons, move the interest of scientists from single neuron 
to small or large networks models. Hence the time seems ripe to summarize 
the contribution of single neuron models to our knowledge of neuronal coding. 
Various types of spiking neuron models exist, with different levels of details 
in the description. They range from biophysical ones on the lines of the clas- 
sic paper of 1952 by Hodgkin and Huxley [SS], to the "integrate and fire" 
variants (see for example [12], [53], [12]). Integrate and Fire (IF) type mod- 
els disregard biological details, that are accounted for through a stochastic 
term, to focus on causal relationships in neuronal dynamics. Their relative 
simplicity make them good candidates for the study of networks. Recent re- 
views discuss qualitative (cf. [H], [TH]) and quantitative (cf. [7D]) features of 
stochastic Leaky Integrate and Fire (LIF) models. These models are a variant 
of IF models where the spontaneous membrane decay is introduced. An older 
paper ([lUZ]) concerns mathematical methods for their study. The aim of this 
work is to collect the existing mathematical methods for LIF models, to pro- 
vide a set of methodologies for future studies on networks. Indeed, although 
the stochastic LIF models are simplified representations of the cells, they are 
considered good descriptors of the neuron spiking activity (cf. for example 
[55] . [73]). Though some criticisms have appeared, showing some lacks in the 
fit of experimental data (cf. [70]), these models are still largely employed. 
The most used is the Ornstein-Uhlenbeck (OU) version but all of them have 
played a role for the understanding of the mechanisms involved in neuronal 
code 

The first IF models date back to 1907, when Lapique ([87]) proposed to 
describe the membrane potential evolution of a neuron, subject to an input, 
using the time derivative of the law for the capacitance. In the presence of 
an input current, the membrane voltage increases until it reaches a constant 
threshold S. Then a spike occurs and the voltage is reset to its resting value, 
to start again to evolve (cf. [134) ). Although it reasonably fitted some ex- 
perimental data this model lied disregarded till the second half of the last 
century. Then it became the embrional idea for "integrate and fire" models. 
The leading idea in the formulation of stochastic IF and LIF models was to 
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partition the features of the neuron in two groups: the first ones were ac- 
counted for by the mathematical description of the neuronal (deterministic) 
dynamics and the second ones globally considered by means of a noise term. 
In Sect. |3]we derive the most popular LIF models after a brief description of 
the biological features of interest in Sect. [5] 

Improvements of LIF models were proposed in the eighties, following the 
initial illusion to become able to recognize the main laws governing our brain. 
The lack of suitable mathematical instruments made soon clear the difficulty 
to determine explicit expressions for the input-output relationship. The end 
of the eighties and the starting years of the nineties are characterized by 
mathematical and numerical advances, accompanied by the development of 
new faster computers. Section |4] is devoted to a review of the main mathe- 
matical methods for the study of stochastic LIF models, updating previous 
reviews [T], [TUS] and [TU7| . 

In the nineties the use of such methodologies, as well as specific reliable 
and powerful simulation methods, allowed to obtain a deeper knowledge of 
the models' features. Unexpected results on the role of noise in neuronal 
coding have been proved mathematically and confirmed experimentally (cf. 
for example |128j ). 

Surprisingly, all research on LIF models has disregarded for a long time 
their ability to fit real data. The only exception was [77] that considered the 
parameter estimation problem. Recently papers on the statistical estimation 
of model parameters started to appear. Section[5]is dedicated to this subject. 



2 Biological features of the neuron 

A comprehensive description of the physiological properties of neurons is out- 
side the aims of this work. We refer to [134] , [135] and [42] for an exhaustive 
exposition of neurobiological properties relevant in the modeling context. 

The neurons are the elementary processing units in the central nervous 
system, interconnected with intricate patterns. Neurons of different sizes and 
shape, but sharing some fundamental features, exist in all the areas of the 
brain. Their estimated number in the human brain is around 10^^. A typical 
neuron can be divided into three distinct parts called dendrites, axon and 
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soma. The dendrites play the role of the input device collecting signals from 
other neurons and transmitting them to the soma. The soma is the non-linear 
processing unit of the neuron. It generates a signal, known as spike or action 
potential, if the total amount of inputs exceeds a certain threshold. The axon 
is the output device carrying the signal to the other neurons. 

The action potentials are electrical pulses, having a duration of about 1 — 2 
ms and an amplitude of around 100 mV. They do not change their shape 
along the transmission. A neuron cannot elicit a second spike immediately 
after a first one has occurred due to the presence of a refractory period. A 
chain of action potentials emitted by a single neuron is called a spike train, a 
series of similar events occurring either at regularly spaced instants of time or 
more randomly. The time between two consecutive spikes is called interspike 
interval (ISI). 

The site where the axon of a presynaptic neuron is linked with the den- 
drite or the soma of a postsynaptic cell is the synapse, which often in the 
vertebrates is of chemical type. When an action potential reaches a synapse, 
it triggers complex bio-chemical reactions leading to the release of a neuro- 
transmitter and the opening of specific ionic channels on the membrane. The 
ion influx leads to a change in the potential value at the postsynaptic site 
and the translation of the chemical signal into an electrical one. This voltage 
response is called the postsynaptic potential (PSP). 

The effect of a spike on the postsynaptic neuron is measured in terms of 
the potential difference between the interior of the cell and its surroundings, 
called membrane potential. In the absence of spike inputs, the cell is at a 
resting level of about —65 mV. If the change in membrane potential is positive 
the synapse is excitatory and induces a negative depolarization, otherwise it 
is inhibitory and hyperpolarizes the cell. In the absence of inputs, i.e. in the 
silent state, the neuron membrane potential decays exponentially toward the 
resting level. 

The dimensions and number of synapses vary for different neurons. Some 
neurons, such as Purkinje cerebellar cells, pyramidal neurons and interneu- 
rons recorded in vitro (cf. for example |70| ) . have a huge number of synapses 
and extended dendritic trees. Integrate and fire models can then be employed 
for the description of their output behavior since, due to the large number of 
synapses, limit theorems can be used (cf. [70], |101j ). 
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3 One dimensional Stochastic Integrate and Fire Models 



3.1 Introduction and Notations 

The huge number of synapses impinging on the neuron determines a stochas- 
ticity in the activating current not considered in the Lapique model. The 
first attempt to formulate a stochastic IF model is due to Gerstein and Man- 
delbrot. In [IT] they fitted a number of recorded ISFs through the Inverse 
Gaussian (IG) distribution, i.e. the first passage time distribution of a Wiener 
process through a constant boundary S. They described the membrane po- 
tential dynamics preceeding the release of a spike through a Wiener process. 
To get a renewal process they assumed that after each spike the membrane 
potential is instantaneously reset to its initial value (cf. for example [27] for 
an introduction on these processes). This model is the basis of successive 
more realistic models. 

In the models, classified as stochastic IF or LIF, one describes the time 
evolution of the membrane potential by means of a suitable stochastic process 
X = {Xt,t > to} with Xto = Xq and identifies the ISI's with the random 
variable (r.v.) first passage time (FPT) of X through the threshold S: 

T^Ts = inf {t>to: Xt> S}. (1) 
The probability density function (pdf) of T, when it exists, is 

g{t)^g{S,t\xo,to)^-^^PiT<t). (2) 

When to = we simply write g{S, t jxo ). In some instances S = S(t). 

In Subsections 3.2, 3.5 and 3.6 we focus on models that describe the sub- 
threshold membrane potential as a diffusion process. In Subsects. 3.3 and 
3.4 we present two continuos time Markov models, the Randomized Random 
Walk and Stein's model. Reviews on IF and LIF models have already ap- 
peared (cf. |105| . |107) . [82) ) but here we unify the notations and we list in a 
single contribution the mathematical results sparse in different papers. 

In the case of models using a diffusion process X — {Xt^t > to}, the 
diffusion interval is / = (^i'')i the drift coefficient and infinitesimal variance 
(the infinitesimal moments) are: 
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with AXt = Xt+At-Xt. The transition pdf / (cc, t \xo, to ) = ^^'^^''^^jf ^ 
is solution of the Kolmogorov equation(cf. [TOT): 



df {x,t\xo,to) df {x,t\xo,to) a'^ixo) d'^f {x,t\xo,to) 

57 +M(2;o) ^ \ ^ TT^ =0 (4) 

oto oxo 2 oxq 

and of the Fokker-Planck equation 
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^^^"^'gf"'^"^ = M^)f {x,t\xo,to)} + \-^ {<j'{x)f {x,t\xo,to)} 

(5) 

with initial delta condition 



lim f {x,t\xo,to) = 6{x - xo). (6) 

to— ft 

Here S denotes the Dirac delta function. 

We suppose that the infinitesimal moments verify some mild conditions 
(cf . [72] , |101| , |107) ) to guarantee the existence of the solutions of the Fokker 
Planck and Kolmogorov equations. Furthermore, when a dependence of the 
diffusion coefficients from t is not specified, the processes are time homoge- 
neous, i.e. their properties are invariant with respect to time shifts. When eq. 
(j4]) is solved in the presence of an absorbing boundary in x = 5, a further 
absorption condition must be imposed: 

\im r{x,t\xo,to) = 0. (7) 

where f°-{x,t\xQ,tQ) ~ -g^P{Xt < x,Ts > t\Xs — y) is the corresponding 
transition pdf. To get renewal processes, X is always reset to xq after each 
spike. 

To characterize a diffusion model, one can also make use of the Ito- 
type stochastic differential equation (SDK) verified by the process (cf. SU- 
SANNE). 



Leaky Integrate and Fire models 



7 



Jump diffusion models, allowing to distinguish the effect of neuronal inputs 
according with their frequency and their size, are presented in Subsection l3.8l 

The role of the threshold shape is illustrated in Subsect. 13. 9[ and the most 
recently introduced IF models are surveyed in Subsection lS.lOl 

To switch from the description of the spike times of the neuron to the count 
of the number of spikes up to a given time t, we introduce, in Subsection l3.11l 
the return processes. 

3.2 Wiener Process Model 

Gerstein and Mandelbrot (cf. described the time evolution of the sub- 
threshold membrane potential through a Wiener process Wt characterized 
by infinitesimal moments 



with e M, : CT > 0. Their model was motivated by experimental observa- 
tions of the ISI's exhibiting histograms typical of stable distributions. Indeed 
this property is exhibited by the FPT of a Wiener process. One gets such 
process from the standard Wiener process W (cf. SUSANNE) through the 
transformation 



To relate the use of the Wiener process with the membrane potential 
evolution, Gerstein and Mandelbrot observed that the Wiener process is the 
continuous limit of a random walk (cf. SUSANNE). The occurring of jumps 
models the incoming of PSP's. The continuos limit is a good approximation 
when the inputs are of small size and frequent . The transition pdf of W is 



/i (x) ~ n (x) = 



(8) 



(9) 



fw {x,t\xo,to) 



dP{Wt<x\Wt„^xo) 



(10) 



dx 
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To mimic the spiking times a constant absorbing boundary S is introduced. 
The spike times are then identified with the FPT, T, of the Wiener process 
originated at Wt^ — xq tlirough the boundary. To obtain the renewal prop- 
erty, the process is instantaneously reset at a;o after each spike. Hence the 
ISI's correspond to the iid r.v.'s r„, n = 1,2, with T„ ~ T. 

The transition pdf of W, if Wto = xo is Gaussian with mean E{Wt) = 
Xq + lJ.t and variance Var{Wt) — a^t, while the FPT pdf through a constant 
boundary S* > xq is an IG distribution, hence the pdf and the crunulative 
distribution are: 



gis,t\xo] 



S - Xq 



exp ■ 



{S ~ Xq- fit) 

2aH 



(11) 



P{T<t)^-{ Erfc 



S — Xq — fit 



2m(S-xo) 



-Erfc 



S — Xq + fJ,t 



aV2t 



(12) 

Here Erfc denotes the complementary error function (cf. [2^). The mean 
value and the variance of the FPT are 



S — Xq (S — Xo)a'^ 

E{T) = -■ Var{T)^^ ' 



(13) 



fi fi" 
The transition pdf in the presence of a constant absorbing boundary S is (cf. 

my- 



/ {x,t\y,s) = . , , exp 

cr^27r(t-s)V 2CT^(t-s) 

exp[-^(5 -y) ^-^^^ ] ) . (14) 



Despite the excellent fitting with some experimental data, Gerstein and 
Mandelbrot model was criticized for its biological simplifications (cf. |135j ). 
However it allows to obtain results that help the intuition for more realistic 
models and it is still used for this aim, taking advantage of the existence of 
a closed form FPT pdf through a constant boundary. 

Its FPT pdf is known also through particular time dependent boundaries. 
These FPT's can be used to account for the refractory period following a 
spike. Indeed a time varying boundary, assuming high values at small times 
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and then decreasing, makes short ISPs rare. The FPT pdf is known for a 
continuos piecewise-hnear boundary S(t) — ai + Pit,: t G i > 1 

where io < < ^2 < ■•• and ai, Pi ^ R with to > 0. If t G [0,oo), such 
boundary is linear and the FPT pdf is 

/ , o . ,\ ^ l«i - xo\ [ai + f3it - fit ~ xpf 
g{ai + l3it,t\xo) = -==r exp . (15 

In the general case, setting a^+i — + Piti one gets that t i— !■ S{t) is 
continuous on [ta, oo). If we put Si = S{ti), the transition pdf for the process 
W in the presence of absorbing boundary S (t), f°'{xi, ti, X2, t2', x, t\xf), to), 
is for t € {tn-i,tn) (cf. [Hg): 



gn 

r{xi,ti;x2,t2; ...■,x,t\xo,to) = ^ ^ — (16) 

{P{Wt, <xi,...,Wt„_, < x„^i,Wt < x;T > t\Wt„ ^ XQ < ai)} 

ri-l 



ri-l 

X 



" ^ ' „(£i-fi)(£i- 

7~. T~ 



n 1 



X 

i=l 



f°-{x,t\Xn-l,tn-l) f {Xi-Xi-i)'' 

X — -^^^^=^^^ exp ' 



V2^(t, - t,-i) V 2(i,-i,_i) 

for Xi < Si, 1 < i < n; xq < Sq and x G (— oo, 5). 

Further closed form expressions for the FPT of a Wiener process have been 
obtained by the method of images (cf. |29l) or as solutions of suitable integral 
equations ( |104j ). when ^^^^ < Ct^", with a < l/2and C a constant. This 
last case is discussed in Section 14.1.1 1 and involves series of multiple integrals. 



3.3 Randomized Random Walk Model 

In the Randomized Random Walk (RRW) the regularly spaced intertimes 
of the random walk between PSPs are substituted with exponentially dis- 
tributed intertimes of parameters A"*" and A~ for excitatory and inhibitory 
PSP's respectively. The process X with Xq — has mean and variance 
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E {Xt) = S {X+ - X-) t; Var{Xt) = S^{X+ + X-)t, (17) 

respectively. Here (5 > is the constant amplitude of PSP's. The FPT pdf 
through the boundary S, with S an integer multiple of 6, is (cf. [135] ): 

g{S,t\0)^ji^—j Is/s{2t^>^y t>0 (18) 

where (.) is the modified Bessel function of parameter rj (cf. f2J). The mean 
and variance of the ISI distribution are (cf. [135j ): 

EiT)^ ^^^J ^ ■ Var{T)^ S {X+ + X-) 
6{X+-X-) 5{X+~X-f 

When (5—7-0, assuming A+ ^ ^ this model converges to the Wiener 
process with /n = 0. 



3.4 Stein's model 

In 1965 Stein (cf. \122[) formulated the first LIF model, i.e. an IF model with 
the leakage feature, by introducing the spontaneous membrane decay in the 
absence of PSP's in the RRW model. The process X is solution of the SDE 

dXt = (^-^ +p^dt + 6+dN+ + S-dN+; Xt„ = xq. (20) 

Here > is the membrane time constant, p is the resting potential, Nj^ and 
are independent Poisson processes of parameters A"*" and A~ respectively 
and (5+ > 0, (5^ < are the amplitudes of excitatory and inhibitory PSP's. 
Generally for this model and for all those descending from it one assumes 
p — 0, since the case p 7^ can be obtained by the shift X 1-^ X — p. Following 
the IF models structure, the spike times are the first crossing times of the 
process through the boundary and the membrane potential is instantaneously 
reset to its resting value after each spike. The infinitesimal moments are: 
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Ml (a:) 



lim 



E[Xt+h~Xt\Xt^x] 



-+ p + \+5+ + \-5- (21) 



h 



E {Xt+h-Xtf\Xt=x 



A+ + {5-f. 



M2 {x) 



lim 

h^O 



h 



The mean trajectory, in the absence of a threshold, is 



E {Xt \Xq = 2:0)= xoe 




) ■ (22) 



The FPT problem for the process pop is still unsolved and the use of simu- 
lation techniques is required for its analysis. 

3.5 Ornstein-Uhlenbeck Diffusion Model 

The OU process was proposed as a continuos limit of the Stein model to 
facilitate the solution of the FPT problem. The rationale for this limit is the 
huge number of synapses characterizing certain neurons such as the Purkinje 
cells. The PSP's determine frequent small jumps and limit theorems can be 
applied to get a diffusion process. The OU process, already known in the 
Physics literature (cf. [136] ). belongs to both classes of Markov and Gaussian 
stochastic processes. 

Different approaches can be followed to obtain the diffusion limit of Stein's 
model. In [71] the convergence of the measure of a Stein's process to that of 
the OU one is studied. In [20], |107j it is proved that the continuos limit of 
the transition pdf for the process (|20|) converges to a pdf that verifies the 
Fokker Plank equation for the OU process. Alternatively the OU model can 
be derived from a differential equation describing the membrane potential 
dynamics. We sketch in the following these last two approaches. 

Due to the time continuity and to the Markov property of Stein's process, 
the Smolukowsky equation holds for the transition pdf: 



f{x,t + At\xo,to)= / f{x,t + At\z,t)f{z,t\xo,to)dz. (23) 



In the absence of inputs to the neuron, the process ([20| . initially in the 
state z at time t, decays exponentially to the zero resting potential, reaching 
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at time t + At the value xi = ze . In the case of an excitatory or 
an inhibitory input at time w G (t, ^ + A£) the potential becomes X2 (u) — 
ze-^t/e_^g+^-(t+At^u)/e ^ ^^-At/e ^g-^-{t+At-u)/e ^ respectively. 

Setting Xi = ^ Jt^^^ (") du, i = 2,3 the l.h.s. of (|23)) becomes: 



f{x,t + At \z,t)=[l- (A+ + A") At] d{x-xi) 

+ \+AtS {x - X2) + X^AtS {x - X3) + o{At). (24) 

Hence ([^ becomes: 

f{x,t + At\xo,to)= e^*/^ { [1 - (A+ + A- ) Z\t] / (xe^*/^ , i |xo , to ) 

+ X+Atf fxe^'/o - ^ (e^*/*^ - l) ,t |a;o,to) (25) 



+ X-Atfixe^ - -^{e^ - l),t\xQ,to) j + o{At). 
Approximating e"^*/^ « ^ + 1, dividing by At, when z4t — > 0, ([26]) becomes 



£?/(2;,i|xo,to) _ d (X 



dt 



— (^/(x,t|xo,to)) +A+ [/(a:;-(5+,t|a;o,to) 



-/ (x,i 1x0,^0 )] + A [f{x-5 ,t\x^,t^) - f {x,t\xQ,tQ)\ . (26) 
Developing the terms in square brackets as Taylor series around x 



df {x,t\xo,tQ 
dt 



d_ 

dx 



[-[^ + X+6+ + X-5-)j{x,t\xo,to) 



(27) 



00 / ^ On 



Assuming (5+ = —6' 
positive constants, for (5 — > we get 



+ A- 



with a,A+,A- 



2(5 



a/ (x,t 1x0,^0) 9 



9f 



dx 



+ m) / (a;,* ko,to) 



2 dx"^ 



(28) 
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i.e. the Fokker Plank equation for an OU process with /i = ^4+ — . Its 
solution with an initial delta condition ([BJ is the transition pdf 



}ou [x,t\xo,to) = 



dx 



2(t-to) 



X exp < 



(t-tg) (*-*q) , 

x — xqc 5 — iia(l — e 5 



2(t^to) 

1 - e 5 — 



(29) 



The diffusion interval coincides with the real line and the mean and vari- 
ance of the OU process Xt with Xtg — xo are: 



Ei 



' {X,\xo) = /.^ (l - e-*/«) + xoe-*/^ Var {X,\xo) = ^ (l " e~''/') . 

(30) 

Properties of the models, as well as the range of validity of some approxi- 
mate formulae for the FPT problem, depend upon the value of the asymptotic 
mean depolarization of the process X. Hence in various instances we distin- 
guish between two distinct firing regimes, subthreshold if E{Xao) < S and 
suprathreshold in the opposite case. 

Model of LIF type can be interpreted in the frame of threshold detectors 
theory. The presence of a feeble noise helps the detection of the signal, a char- 
acteristic of any threshold detector. In Fig. [T] we plot the mean value of an OU 
process together with E{Xt) - War{Xt) and E{Xt) + 3Var{Xt), mak- 

ing use of (pOb ) . The two panels correspond to examples in the subthreshold 
regime (Panel A) and in the suprathreshold regime (Panel B). The intrinsic 
random variability determines crossings even in the subthreshold regime. 

The OU model can also be obtained from the differential equation for the 
time evolution of the subthreshold membrane potential in the presence of 
spontaneous decay of parameter and net input /x: 

dXt Xt , , 

— ^--+^,X,^x,. (31) 

Adding a noise term of intensity a to account for the random PSP's, one 
gets: 
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Fig. 1 Mean value (middle line) and curves E(Xt) — ZVar(Xt) {lower line) and E(Xt) + 
'iVar{Xt) {upper line) for an OU process with parameters ^ = 0.8 mVms~^ , = 0.2 
mV^ms~^ and 9 = 10 ms. 



dXt = ( +^JL]dt + adWt; Xq = xq. 



(32) 



This is the SDE of an OU process (cf. SUSANNE). The analytical expres- 
sion for the FPT pdf of the OU process is still an open problem. In [5] three 
alternative representations of the distribution of T are introduced. The first 
one involves an eigenvalue expansion in terms of zeros of the parabolic cylin- 
der functions while the second is an integral representation involving special 
functions. We report here the third one that writes the FPT pdf of an OU 
process with fi — through the boundary S in terms of a three-dimensional 
Bessel bridge: 



lit) 



gw{t)EBb 



exp 



/ (rs-Srds 



(33) 



Here gw{t) is the FPT pdf through the boundary S for the standard Wiener 
process, is the three-dimensional Bessel bridge over the interval [0,t] be- 
tween To = and rt — S ~ xq. This process is solution of: 



dr. = 



y - rs 
t-s 



1 



ds + dWs, tq = X, s < t. 



(34) 
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In (|33p Esb indicates the expectation with respect to the Bessel bridge law. 
In [138| formula is used to approximate the FPT pdf with Monte Carlo 
techniques. 

An explicit expression for the FPT density of continuous Gaussian pro- 
cesses to a general boundary is obtained under mild conditions in [37_ , while 
an expression for the FPT of the Wiener process to a curved boundary is 
expanded as a series of multiple integrals in [35] . 

Existing available closed form expression include the case of a hyperbolic 
boundary (cf. [Tl] ) 



with A and B arbitrary constants (cf. |101| V Furthermore specific boundaries 
can be obtained through the space time transformations described in Section 
14. 1.21 applied to closed form solutions for the case of the Wiener process. The 
Laplace transform of the FPT pdf in the case of a constant boundary S is 



where D^, (.) is the Parabolic Cylinder Function (cf. [2 ) of parameter i/. No 
analytical inversion formula is available for eq. (p6|. Reliable and efhcient 
procedures, discussed in Sect. |4l can be applied to obtain the FPT pdf ei- 
ther numerically or by means of simulations for constant or time dependent 
boundaries. 

The FPT mean has been determined as derivative of (pS)) . computed in 



(35) 



(cf. m) ■■ 




(36) 



A = (cf. [T02] l: 




(37) 
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where xi = {^9 — Xo)\/2/ [a'^O), xs ~ {fJ-0 — 5)-\/2/ (a^O) and (j){a,c]z) is 
the Kummer function (cf. ^2,). Ahernatively the mean is expressed through 
the Siegert formula (cf. |130| ): 

fire f^"''^ r z 1 

EiT) . ^- J_^^ {l + ^^/(^)|exp(^)^., (38) 

where Erf{.) denotes the error function (cf. [5]). 

Use of ((37l) or ((38)) depends on the value of the parameters since the 
two formulae present numerical difficulties for different ranges. Approximate 
formulae (cf. 84 ) hold for specific ranges, li jiO > S and cr — >■ 0, i.e. in the 
quasi-deterministic case, the mean FPT can be approximated by equating 
the expression of E{Xt\xQ) with S to obtain (cf. 81 ): 

^m^-oK^)- '''' 

Note that (p9| disregards the effect of the noise on the crossings. If a;o << S, 
or equivalently if a is sufficiently small and /i is negative so that the crossing 
is rare event, the approximation 

holds (cf. |47j). A linear approximation for the firing rate / 
tained using is (cf. [HI]): 

f{p) = ^(cjV^ +29^^-3). (41) 

This approximation holds when /a and (^^\/9 — S9^ /a are small 

enough. 

When neither (|37|) nor (|38l) are suitable for computations and /ir ^ S* but 
a is not small enough to apply approximation (j39p an "ad hoc" procedure 
to evaluate the mean FPT is possible. One establishes at first the time ti at 
which E{Xt) + 2.yVar{Xt) crosses the threshold 5, i.e. when most trajecto- 
ries are still below the threshold. For t > ti the process is then approximated 
by means of the Wiener process with drift /x and initial value E{Xti). 

The second moment of the FPT for the OU process is (cf. [92j): 



(40) 

= l/E{T), oh- 
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E{T'^) = 2eE{T) 

+ 26*2 I 7^ In 2 



where xi, xs are defined as in ((57)) and 
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(42) 



^2fe+l 



fc=0 
2n+3 



^^(n+l)!(2n + 3)^^ 



/c!(2fc + l)' 

^ 2fc + 1 



V'i(z) = 2 / e- / e 
/o Jo 



2fcz2'=+2 



fc=0 



(2A: + l)!!(fc+l)' 



V'2(z) = J2 



2n^2ri+4 

^0 (2-T3)!!(^+ 2) 



" 1 

fc + 1 



(43) 



In [22] the mean, variance and skewness of the FPT for the OU process 
are tabulated for neurobiologically compatible choices of the parameters. 
Asymptotic results for the FPT of the OU process are presented in Secion 



3.6 Reversal Potential Models 

The diffusion interval of the OU process is the real line but large negative 
values of the membrane potential are unrealistic. Hence other models intro- 
duce a saturation effect on the membrane sensibility. When the value of the 
membrane potential is close to the reversal potential Vj the incoming inputs 
produce a reduced effect (cf. [75], [M]). A diffusion model with reversal po- 
tentials is proposed in [44] as a diffusion limit on a birth and death process. 
A similar diffusion limit is obtained in [78 from a variant of Stein's model 
(|20p where an inhibitory reversal potential Vj is introduced: 
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1, 



dYf = --Yfdt + 



e {Yt - Vi) + iyjYt - Vi dN,-; Yo = yo- (44) 



Here N^, NjT and 6 are the same as in (PO)) . e € (—1,0), Vi < xq are 
two constants and ^ is a suitably defined random variable. The first two 
infinitesimal moments ([3|) of this model are: 



fi{y) = -y(^^-eX-j+S+X+-eX-Vi; (45) 
a^iy) - A+ (5+)' + X-e^ (y - Vjf + X-Var (0 {y - Vj) . 
The mean trajectory of the process originated in Yq — yo is 

EiY \yo) = xoe-(i--) + 'T'f (l - e-(i--)) . (46) 

The diffusion limit of this model (known in the neurobiological literature 
as the Feller model and in other contexts as the Cox-Ingersol-Ross process) 
is identified with the solution to the SDE (cf. [78] ) 

dYt - f - — + M2 ) dt + G2^JYt-VldWu Y^^y^. (47) 



T 

Here the constants [ii , a2 and r are related with those of model (|44)) by 

imposing the equality of the infinitesimal moments. One first sets r — 

and /i2 = A+(5+ — eX^Vj = A+(5+ — X^ where Aj — eVi, then substitutes 

the variable ^ in (jH]) by a suitable sequence of r.v.'s such that in the limit 

for n — > oo one gets Var{£^n) = 0. This choice allows to obtain the same 

infinitesimal variance at the resting level for the two processes. Hence one 

gets 

^2 = • (48) 

Note that, due to the expressions for r and /X2, the parameters t and 112 
appearing in the SDE (|47|) bear a different meaning here with respect to the 
corresponding parameters 9 and ^ in the OU model. Furthermore t < 9. 

The diffusion coefficient of the process (l47t becomes negative if Xt < Vj, 
hence the diffusion interval is / = [V/,00). The boundary in Vj is regular 
or exit, depeding upon the values of fi2 and (72, according with the Feller 
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classification of boundaries (cf . 72 ) . To determine tlie transition probability 
density of the process (H71) a boundary condition should be added. The natural 
choice, that respects the model features, is the reflecting condition: 

\im^ l^fi{x)f {x,t\xo,h) - -^[a^{x)f {x,t\xo,to)]^ ^0. (49) 
The Feller process is generally known in its standardized form 

dXt =(-—+ /iF ) dt + a^^fXtdWt; Xo = -Vi, (50) 



T 

that can be easily obtained from (I47|) by performing the space transformation 
— > — Vi and by setting fip = /i2 — This equation is defined over 
/ = [0,oo). 

A further notation for the parameters of the Feller process, largely em- 
ployed in the literature, sets: 

1 '^i 
p = ~-; q = ^ip; r= —. (51) 

T Z 

The transition pdf of the Feller process X depends (cf. [77) upon the 
nature of the lower boundary in x — for the process (|50p or in x = Vj for 
the process (|47p and on the selected boundary condition for the solution of its 
Kolmogorov equation. If we impose a zero-flux condition ()49p in the origin, 
using the notation (|5ip . we obtain the transition pdf (cf. [44]): 



fpe ix,t\xo,to) 



dPiXt<x\Xt„^xo 



dx 



X cxp ■ 



r(eP(*~*o) — 1) 
r (eP(*-*o) - l) J r (eP(*^*o) - l) 



(52) 



Here In{z) indicates the modified Bessel function of the first kind (cf. ^) of 
parameter r]. 

The mean trajectory of the process (HTl) originated in Xq = xq is 



E {Xt\xo) = xae-'l^ + ^i2T (l - e-'/") 
while its variance is 



(53) 



20 Laura Sacerdote and Maria Teresa Giraudo 

Var{Xt\xo)^Tal{l-e-i) | ^'^~^' (l - e^) + (xo - e'^ | . 

(54) 

The FPT pdf of the Feller process cannot be obtained in a closed form 
but it can be evaluated by employing the methods described in Sect. S) Fur- 
thermore its Laplace transform is (cf.|43j): 

(p ( ^ 3.- P^o \ 
gx{S\xo)= )\"" (55) 

\ p ^ r ^ T j 

Here denotes the Kummer function (cf. [2]) and the notation in (|5ip has 
been employed. The mean firing time, when a;o < S" < Ve, is (cf. [S]): 



A^2r - V, \- ' ^^{n + 1) nLi (/^2r -Vi+ zra|/2) ) ' 

(56) 

If r/i2 » S and CT2 is suitably small, the mean FPT can be approximated 
with a formula analoguous to ((5^ for the OU process. Indeed it holds: 

EiT)^-rln(^^]. (57) 



When the crossing is a rare event, i.e. xq « S or (T2 is small, a result 
analogous to (^(11) can be derived (cf. |47)): 



2(m2T-V-j) 

- 4 



L2(5-y/) 



(58) 



Here denotes the Gamma function (cf. [2]). 

The second moment of the FPT has been obtained in 
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2tE{T) 



s-vj+y: 



(S-Vi) 



i(fc+i)n.:=i 



2t^ [{S-Vj f+'-ixo~Vj)^+'] 

itl2T - Vl)ik + 1) 




^J■2r - Vi 



k+l 



(59) 



The expression for the third moment can be found in |107) . 
A further variant of Stein's model was proposed in [SD] with two reversal 
potentials, an inhibitory (lower) one Vi and an excitatory (upper) one Ve'- 



dXt = -^Xtdt + S+{Ve - Xt)dN+ + [5-{Xt - Vi) 
- J^{VE-Xt){Xt~Vi)\ dN[- 

Xo = Xq. 



(60) 



Here the two independent Poisson processes and have intensities 
A"*" and A~, respectively, —1<5^ <0<(5+<l,Jisa r.v. defined over 
the interval (—1 — 6^,-6^) such that E{J) = 0. For a sequence of models 
(|60p indexed by n one can assume that — > 0- , , — > cx) in such 
a way that (5+A+ — > ^ > 0, (5,7 A^ v < 0. Simultaneously E{J^^) — > 0+ 
in such a way that X^E{J^J 
approximation 



cr| > 0. This allows to obtain the diffusion 



dXt = (-^ + Ms) dt + <j3^{Xt-Vi){VE-Xt)dWu Xo = ^0 (61) 

where the new constants T3 and /is — (xVe — vYi have been introduced. 

The diffusion interval of this process is [V7,yE] and its transition pdf is 
solution of the Fokker-Planck equation with initial delta condition and with 
reflecting conditions of the type (H^ on both boundaries Vi and Ve- 

The first two moments of the membrane potential (cf. 80!) are: 
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E{Xt \xo) = xoe-*/^^ + M3T3(1 - e-*'/^') 



/3(2/3 + a|) _„,(/3~a)(2/3 + a2) 



K2a + a|) 



a(a + cr|) 



+e 

[213 + a' 



(2a + f7|)(a + t7|) 
3 , „_.t_.,^*(2a-2/3 + a|) 



1 - 



xo - Vi 



Ve-Vi 



a + cr.3 





' 2:0 


-VI^ 






-Vi\ 



Here a = I/T3 and /3 = (— aV/ + ^3) / (Ve — Vi) and a typo in [5D] is cor- 
rected. Use of (|62|) and (|63)) allows the computation of yar(Xt |a;o ). 
The mean firing time through a boundary S < Ve is: 



E{T) 



S - Xq 
I3{Ve - Vi) 

^ r{2Plal + l)r(2a/a| + 1) (5 - V,)"+^ - (rrp - 
r(2/3/a| + n + 2)r(2a/ai) /? (n + 2) (1/^, - 



(64) 



If the boundary crossing is a rare event, a result analoguous to (|^ and 
Ml) holds ([80]): 



i?(r) 



/3(V^B-'^^/) 



a{S + xo-2Vi) 
{2p + al){VE-Vi) 



When ^3X3 > 5, one gets a result analogous to and ([57| : 

S - T3^l3 



E{T) « -T3 In 



2^0 - T3M3 



(65) 



(66) 



3.7 Comparison between different LIF models 

The mathematical complexity of the FPT problem increases with the at- 
tempts to make the models more realistic. However it is desiderable to avoid 
the use of complex models when they do not add any improvement with 
respect to the simpler ones. 

In Fig. [5] we compare sample paths from the OU, the Feller process and the 
process with double reversal potential, simulated using Euler algorithm (cf. 
SUSANNE) on the same leading Wiener process trajectory. Furthermore 
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we consider 6 ~ t ^ 10 ms and we choose the same level of variability at the 
resting level for all models. Hence cr| = cr^ and cr| Vg; | V/ 1 = cr^ . Since the 
three processes do not show relevant discrepancies, this analysis suggests to 
prefer the first two models due to their better computational tractability. 



30 



~ 2H 40 60 BO 100 

t (ms) 

Fig. 2 Sample paths of the OU {dashed line), the Feller {dotted line) and the double 
reversal potential (continuous line) models employing the same leading Wiener process 
realization. Here ^ = /i2 = /^s = 1 mVms~^ , = 0.9 mV^ms~^ , Vj = —10 mV, 
Ve = 30 mV. 

When one wish to compare the FPT pdf's one gets different results ac- 
cording with the selected criterium for the parameters values. In [81] the OU 
and of the Feller ISIs, computed through (??), are compared, according with 
three different criteria: 

• to get the same values for their corresponding discrete versions: 

= M2; CT^ = -(T2 V/, 

9 and r chosen accordingly (cf. Fig. IJK); 

• to get the same mean trajectory for both models (Fig. [2P): 



61 = r; /i = /X2 
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• to get almost equal FPT densities. For this aim, one fixes the parameters 
for one model and determines a set of parameters, reproducing a similar 
ISI distribution, for the other one (of. Fig. [3]C). To guess possible set of 
parameters for the second process we impose the equality of the mean and 
the variance of their FPT's. 

The last case illustrates the evenience where a histogram of experimentally 
obtained ISPs can be fitted by either of the two model distributions. 



A. B. C. 




t (ms) t (ms) t (ms) 

Fig. 3 Comparisons between the OU {continuous line) and the Feller {dashed lines) mod- 
els. Panel A: Vi = —10 mV, 5 = 10 mV, 9 = 6.2 ms; parameters controlling the PSP sizes 
and the intensities of the input processes: = —ij = 2 mV, e = —0.2, A = 8/6 = 1.290 
ms-\ u) = 4/9 ^ 0.641 ms~^. Panel B: xq = mV, Vj = -10 mV, S = 10 mV, 
9 = T = 6.2 mV, cr^ = 4, 9 and 16 mV'^ms~^ (from bottom to top), fj, = fip = 
niVms~^ . Panel C: Feller model yo = mV, 5 = 5 mV, Vj = —10 mV, t = 6.2 ms, 
= mVms~^ , (Tj = 4 mVms~^ ; OU model xq = mV, S = 5 mV , 9 = 6.2 ms, 
= -0.799 mVms-^, ct^ = 43 93 mV^ms''^. 

The use of the more complex Feller model seems preferable when one 
has data of the membrane potential between consecutive spikes ([SI])- When 
only the ISI distribution is available, both models fit the data. In (jll9|) 
qualitative comparisons between the OU and the Feller processes, obtained 
through the stochastic ordering techniques (cf. |123| ) in Subsubsect. 14.1.51 are 
presented. In ( [118] ) the same techniques are used for a sensitivity analysis on 
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the parameters of the FPT pdf. Stochastic ordering properties of the FPT's 
are used in [SD] to select the model. 

Membrane potential data analyzed in |67| show that the same neuron, 
under different experimental conditions, can be described either by the OU 
model, by the Feller model or by an alternative model with a quadratic 
diffusion coefficient. 



3.8 Jump diffusion models 



In Subsections (13.51) and p.6p . to perform the diffusion limits, it was assumed 
that all the contributions to the changes in the membrane potential were of 
the same small amplitude and the frequency was large enough (cf. ;20J). 
However PSP's impinging on the soma can play a different role with respect 
to the contributions on different areas of the neuron. 

LIF models where the subthreshold membrane potential dynamics is de- 
scribed by jump diffusion processes allow to separate inputs according to 
their strongness. Jump diffusion models can be obtained from a variant of 
the Stein-type model: 



j=i k=i 

Xtg = XQ. (67) 

Here N^,Nl are independent Poisson processes of parameters and A* and 
amplitude 6'^ and 6^ accounting for the strong contributions. N^'-' , NfT'^ 
are independent Poisson processes of parameters Xj' and , independent 
from and Nl . H Sj' , ^ and at the same time A^ , A^ — ?► cx) so that 
Sj'X^ + ^J^J ~^ M ^-nd (^j^)^A^ + (^7 -'^'^7 ~^ ^ diffusion approximation 
can be performed to get a process solution of the SDK: 

dXt = ( +fi)dt + a^dWt + S'^dNI + S^dNl; Xt„ = xq (68) 



where W is independent from and N'l . The model ([68l) is a jump diffusion 
process with an OU underlying diffusion. Other jump diffusion models may 
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be obtained introducing the reversal potentials. All these models are of LIF 
type, requiring the superposition of a boundary S to mimic the spike activity. 
The crossings can occur either for diffusion or for an upward jump when 
X £ {S ~ 6'^, S). Hence the spike time is the time of first entrance (FET) into 
the strip {S, oo). The cases of underlying Wiener with drift and OU process 
have been considered in ([79], [S^, |120j . |131j ). The exponential distribution 
for the jump epochs preserves the Markov property of the process In 
|120j and |124| the case of IG distributed jump epochs is discussed but in this 
instance the resulting process is no more a Markov one. 

To compute the ISI distribution for IG and exponential time distributed 
jumps one resorts to simulation techniques. Differently from the unimodal 
behavior of the ISI distribution of diffusion models, jump diffusion ones have 
a multimodal shape (cf. Fig.H]). 

The only analytical results on the FET problem for jump diffusions refer 
to an underlying Wiener process with constant boundary. Lower bounds are 
proposed in [31] for the FET density and in [57] for the FET mean and 
variance, together with exact formulae for the specific case of large jumps, 
when the jumps are driven by a counting process. An approximate solution 
to an integral equation for the FET density of a jump diffusion process is 
discussed in [SI] for the Wiener process. 

3.9 Boundary shapes 

Gonstant thresholds are typically employed in LIF models for their easier 
mathematical tractability. However the refractory period following each spike 
has been modelled by means of threshold shapes (cf. |101j . [52]). These bound- 
aries attain very high values after a spike, then decrease under the reference 
value and finally oscillate around a constant value (cf. Fig. jS] panel A). 

In [23] a dynamic threshold obeying to a differential equation is considered 
for the same aim. A boundary which is a linear combination of two exponen- 
tials with different time constants is proposed in [75], to fit experimental 
data. The use of this boundary, together with the lack of the resetting of 
the membrane potential after a spike, allows a very goood fit of the data. 
A computational method that can reproduce and predict a variety of spike 
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Fig. 4 Examples of ISI distributions for jump-diffusion processes II68I I; underlying OU 
diffusion process with different parameters for Panels A and B, C and D. Panel A: IG 
time distributed jumps. Panel C: Exponential jump interarrivals. Panel B, Panel D: 
ISI distribution for pure OU diffusion, parameters as in Panels A and C respectively. 

responses has been deviced in [76] using a multi-timescale adaptive threshold 
predictor and a nonresetting leaky integrator. 

In [24] , [25] thresholds with fatigue are proposed to account for experimen- 
tal data showing a progressive decrease of excitability during high frequency 
firing. This type of threshold destroys the renewal and Markov character of 
the process but allows to describe adaptation phenomena through LIF mod- 
els. The inclusion of time dependent boundaries prevents the use of many 
mathematical methods described in the next Section, however reliable nu- 
merical and simulation techniques can be applied (cf. Sect.Hl). 

Finally, the study of periodic boundaries or of noisy thresholds ([HS]) be- 
comes a useful mathematical method to deal with periodic inputs. Indeed 
one transforms the original process with time periodic drift and constant 
boundary into a time homogeneous diffusion process, constrained by a pe- 
riodic absorbing boundary (cf. |100j . [129] . [51]). To illustrate this idea let 
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US consider an OU model with periodic input of frequency w, phase Lp and 
amphtude A. The SDE ^ for this LIF model is 

dXt = + Asin{ut + (^)j dt + adWu Xq = xo (69) 

with xq < S. X is not an OU process, however the change of space variable 
AO 



Yt^Xt- 



sin(ajt + ip - ^) — e su\{if) — ^) 



(70) 



with ^ — arctan(aj(?) transforms (|69|) into the SDE of an OU process with 
parameters 9, jj. and ct, t/o = xq, and the constant boundary S into 



S'(t) = S- 



AO 



sin(a;i + (/? - ~ e"*/** sin(v3 - ^) 



(71) 



The ISFs of the periodically modulated LIF model with constant threshold 
are distributed as the ISFs of a LIF model with constant input and appro- 
priate time-dependent threshold (cf. Fig. [5]). 
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Fig. 5 Panel A: Time dependent threshold S{t) = 10 + 3e~* -0.6e Too sin( j^). Panel 
B: Threshold of the transformed process (upper) and simulated FPT pdf (lower) for an 
OU process with additional sinusoidal term in the drift coefficient. 



Leaky Integrate and Fire models 

3.10 Further models 



29 



New efforts on LIF models are mainly devoted to the study of input-output 
relationships or to the analysis of small neuronal networks with units de- 
scribed by LIF models (cf. |131j ). New variants of LIF models have recently 
appeared in the literature to catch further features such as plasticity (cf. [55]) 
or to improve their flexibility and their predictive power (cf. for example |13j , 
[21], [H]). We quote also the LIF compartmental models (cf. [53], [TOH], [TTU] . 

that extend the one dimensional ones by introducing systems of SDK's 
to describe the dynamics of different components of the neuron such as the 
dendritic and the soma zone in the case of two compartmental models. These 
models require the study of the FPT of one component through a boundary 
to describe the ISI's. Up to now this problem can be handled only through 
simulations. 

3.11 Refractoriness and return process models 

An alternative approach to the study of spike trains focus on the number 
of spikes in a prescribed time interval. This study allows to introduce the 
refractoriness of the neuron in a quite natural way. To this aim one can 
associate a return process {Z — Zt,: t > 0} \tv the interval (/, S), S € I, to any 
regular diffusion process X — {Xt : t > 0} on I = {I, r) as follows. Starting 
at xq € {l,S) at time 0, the process Z coincides with X until it attains the 
level S. At this time it is blocked on the boundary S for a random time 
and no new crossings can occur during this refractory period. Then Z and 
X are instantaneously reset at xq and the process Z evolves as X until the 
boundary S is reached again, and so on. We show in Fig. (BJa sample path of 
this process when the refractory period is constant. 

Let Fi and Ri, i = 0,1, be the r.v.'s describing the time between the 
i'th reset and the {i + l)-th crossing and the i-th refractory period. For time- 
homogeneous diffusions, the r.v.'s Fi are iid with pdf g{S,t\xa). It is also 
assumed that the r.v.'s Ri, i = 0, 1, are iid with pdf ip{t) depending only 
on the duration of the refractory period. 
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Fig. 6 Sample path of a return process Z with a constant refractory period. 



A counting process M = \Mt^ i > 0} can be introduced to describe the 
number of attainments of the level S by the process Z up to time t. Let 



|xo ) = P {Mt = fc |Zo = xo } , A: = 0, 1, ... 
be the probability that fc firings occur up to i. Then (cf. [IDS') 



(72) 



go(i |a;o) = 1 - / 5(S',T|xo)dT 
Jo 



g{S,T\xo)dT 



Lp{T)dT 



(73) 

g{S,t\xo) 
fc = l,2,... (74) 



where * means convolution and exponent (fc) denotes (fc)-fold convolution. 

Expressions for such probabilities have been obtained in [Sj for the Wiener 
process for exponentially distributed refractory periods. In the general time 
homogeneous case 



E{Ml\xo} 



k>l 



k'-qk{t\xo) r = 1,2,... 



(75) 



is the r-th order moment of M.Let 7^, i = 0, 1,2, ... be the r.v.'s describing 
the ISFs and let Iq be the time of the first firing. One has (cf. [49 ): 



Leaky Integrate and Fire models 



31 



E{I) = ti{S \xo) + E{R); E{I^) = t2{S l^o) + E{R^) + 2E{R)ti{S \xo ) 



E{I^) = tsiS \xo ) + EiR"^) + 3EiR)t2{S \xo ) + 3E{R^)hiS \xo ] 



(76) 



where tr{S \xq) is the r-th order moment of the FPT. If the first three mo- 
ments of the refractory period are finite, the following asymptotic expressions 
for large times hold for the first two moments of M (cf. [49]): 



E 



E{Mt\xo}^ 
{[Mtf\xo}^ 
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m 

1 E{P) 
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ti{S \xo). 



_ tijS \xo 

m 



(77) 
(78) 



E^I) E^{I) 

In [3S] and |106) the case of absence of refractory period and that of a 
random distribution for the return value are discussed for the Wiener, OU and 
Feller processes. Alternatively [17 proposes to model refractoriness through 
return processes characterized by an elastic boundary as firing threshold. 



4 Mathematical methods for First Passage time problem 
and their application to the study of neuronal models 

We update here previous reviews ([T], |105) and |107j ) on the methods avail- 
able up to now to deal with the FPT problem for stochastic LIF models. 

In the case of diffusion processes, closed form expressions for the transition 
pdf's are determined as solutions of the Kolmogorov or the Fokker-Plank 
equations of Subsections l3.2l[?3] and [?751 The Fourier transform is the typical 
method to get these solutions. 

Closed form solutions for the FPT pdf refers generally to specific time 
dependent boundaries. The constant boundary case is known for the Wiener 
process (cf. ^) or for the OU one if 5' = 0, a* = and 7^ 0. The FPT 
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pdf has been determined for the hyperbohc shape boundary (j35|) for the OU 
process (cf. [lOlJ ) or for the Feller process (cf. |ll5 ,) in the case of bound- 
aries corresponding to symmetries for these processes. Use of the reflection 
principle allows to determine the FPT pdf through Daniels' boundary f[28j. 
US]) [SE]) while a particular FPT pdf is found through a symmetry based 
constructive approach in [45] . The theory of group transformations is used in 
|115] and [53| to determine the transition pdf's of a Feller process between 
the origin and a hyperbolic-type boundary and of an OU process between 
a lower reflecting and an upper absorbing constant boundary, but is of no 
interest for neuronal application and hence we omit the description of this 
method. 

In Subsect. l4.1l we review the commonly used analytical techniques for FPT 
problems: integral equations (14.1.11) . change of variables or of measure (|4.1.2I) . 
asymptotic studies (j4.1.3p . computation of FPT moments (|4.1.4p . stochastic 
ordering (|4.1.5p . In Subsubsect. 14.1.61 we present the available methods for 
jump diffusion processes. In Subsect. (|4.2[) we introduce the direct (|4.2.1I) and 
inverse (|4.2.2p FPT problem. Finally in (14. 3p we sketch specific simulation 
techniques for FPT's and the numerical tools for their solution. 

4-1 Analytical methods 
4.1.1 Integral equations 

In 1943 Fortet ([40]) proved, under mild conditions for the boundary S{t), 
that the Volterra integral equation of the first kind: 



holding for x > S{t), holds also for x — S{t). 

When the boundary is constant and the process is time homogeneous, eq. 
([75]) is of convolution type and the Laplace transform method can be applied. 
Denoting as f\{S \xo ) and gx{S \xo ) the Laplace transforms of f{S, t \xo ) and 
of g[S, t\xQ), for x > S > xq one gets: 




(79) 
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gx{S\xQ) = . 80) 

Generally the Laplace transform cannot be analytically inverted due to its 
complex expression (cf. for example (1361) ). 

Eq. ((79|) for x ~ S(t) has a weakly singular kernel for t ~> t. In- 
deed, any diffusion behaves as a Wiener process for small times. Hence 
/ {S{t), t\S (t) ,t ) w ^^^= with k{t, r) -H- as r — > This makes numerical 
methods for their solution unstable. Hence it is convenient to switch to a 
second type Volterra equation. Integrating (1751) between the left side of the 
diffusion interval I and S{t) and then differentiating with respect to time, one 
gets a second kind Volterra equation (cf. |104| ): 

g{S{t),t \xo ) = 2j{S{t),t \xo)^2 fg {S{t), t \xo ) j {S (t) , t \S (r) , r ) dr, 

(81) 

Here we have introduced the probability current through z at time u of the 
diffusion process X whose pdf is solution of ([5]): 



1 ( d 

j{z,u\w,v) = fi{z)f{z,u\w,v)- - |— [cr^ {y) f {y,u\w,v)] 



■ (82) 



Eq. 1ST] has a weakly singular kernel. For the Wiener process, when 
< Ct'^", with a < 1/2 , C a constant and limt_>o <S'(t) > Wt^ = xq, 
g {S{t), t \xq ) is the only L? solution of (|8T|) . It can be expressed as (cf. [104] ) 



dS{t) 
dt 



g (Sit), t\xo)^ 2j {S{t), t\xo)-A f drj iS{t), t \S{t),t)j {S{t), t \xo ) 



^4" / drjn{S{t),t\S{T),T) 



(83) 

X ^2j(5(t),t|x-o)-4 / d0j{S{r),r\S{0),e)j{S{e),e\xo) 







where 



Jn{S{t),t\S{l 



d0n {S{9),e \s{t),t ) j„_i is{t), t \s{e),e) (84) 
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for n = 2, 3, ... and 



n{S{t),t\SiT),T)= J dejiSie),9\S{T),T)jiS{t),t\Si0),9). (85) 
A third integral equation was proposed in |14) : 



g iS{t),t\xo) ^ -2i^iSit),t\xo) + 2 / g {S{t),t\xo) {S{t),t\S (t) ,t) dr 

Jo 

(86) 

where 



(Sit), t\x,T)^-{F {Sit),t \x, T)} + k{t)f {S{t), t\x,T). (87) 

S(t) 

Here F {S (t) ,t\x,T) — f {y,t\x,T) dy and k [t) is an arbitrary contin- 
uous function. A suitable choice for k (t) allows to make the kernel of (|86p 
regular and hence eq. (|86p becomes optimal for numerical integration. For ex- 
ample the expressions (|87|) for the OU and the Feller processes, respectively, 
that make the kernel of (j86| regular are (cf. [T4] ) 



■)pou (S (t),t\x,T)) = 



.s'{t) + s{t)/e- n 



e e 



{[S{t) ~ ^,e]e 



2(t-T) , 

(1-e^) 



x + tie]]f{S[t),t\x,T)) 



and 



ll^Fel (S (t),t\x,T)) 



S{t) 



r[ePit-^) - 1] 
S'{t) 



■ exp 



p[S{t) + xeP^*-^'>] 



pS'(t)eP(*-^) 1. r 



X /, 



q/r-l 



2py^S{t)xeP<-^-' 
r[ePit—^) - 1] 



p^S{t)xeP(t-) 



2py/S{t)xeP(t-' 



(89) 
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Other choices for k(t) make the integral on the r.h.s. of ([55]) vanish for bound- 
aries with particular symmetry properties, such as psp for the OU process. 
Infinite sum expansions, bounds and approximations for g {S(t), t \xq ) can be 
obtained from ^ (cf. [TT7] ). 

Expressions that regularize the kernel of eq. ([55]) can be found also in other 
specific cases. 



4.1.2 Change of variables or change of measures 

The transition pdf and the FPT pdf of an assigned processes can be obtained 
through changes of variables and/or changes of measure. 

Let Y — {Ft, t > 0} be a diffusion process with / C SR, characterized by 
drift fJ.{y,t) and diffusion coefficient fj{y,t). One may wish to transform this 
process into a Wiener process through suitable space-time transformations, 
when this transformation exists. In [lOOj it is shown that a transformation, 
conserving the probability mass, mapping the Kolmogorov equation of the 
process Y into the analogous equation for the Wiener process 

df d^f 

with initial delta condition is of the form 

r' = 0(T); y' = V(y,T); /(x, t |y, T)da; = /'(x', r' (91) 
This transformation exists if and only if the infinitesimal moments verify 



(92) 

Here z g / is an arbitrary value and ci{t) and C2{t) are arbitrary functions 
of time. If (l92l) holds the transformation is: 
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duci (u) exp 



r' = (/)(t) ~ ki duexp — I 

J Ti L J To 

fix,t\y,T)dx^ f'{x',t'\y\T')dx' 



duc2{u) 

d0c2{e) 

dec2{0) 



y dx 



cr(x,T) 
-fc2 



(93) 



where z ^ I, Ti £ [0,oo) and fc^, i = 1,2,3 are constants with ki > 0. For 
example the transformation 



y' = i^{y,r) 



ki - n) 
(2/ - 2) H 



{r) 



kl9 r 2(t~to) 2(ti-to) 

e 5 — e 5 



f{x,t\y,T)dx = f'{x',t' \y',T')dx' 



(94) 



changes the OU process into a Wiener process. Here ti,T2 > are arbitrary 
times. The transformation ([M)) sends the hnear boundary S{t) = a + bt for 
the Wiener process into the U-shaped boundary ([55]) for the OU process. 

The Feher process can be transformed into the Wiener one, conserving the 
probabihty mass, only if ^ = i. In [51] a necessary and sufficient condition 
analogous to is given to transform a diffusion process Y into a Feller 
process. 

A change between the measures of two processes is considered in [116j , [5] 
and [5] . In [S] the Girsanov theorem (cf. |112) ) and the change of measure 



dP' 



OU 



exp 



1 



Xq 



1 

2^ 



dP 



(95) 



are applied to the OU process to obtain its FPT pdf through a constant 
boundary. Here P'~'^ and P denote the distributions of an OU process with 
11 = and of a standard Wiener process, respectively. 
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4.1.3 Asymptotic results 
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Asymptotic results play an important role in the study of the FPT pdf be- 
cause they hold from relatively small values of the involved variable. Studies 
on the asymptotic behavior of the FPT pdf belong to two different classes: 
large values of the boundary and large times. Let us first list asymptotic re- 
sults in the first class. In (fSl]) the asymptotic exponentiality of the FPT for 
an OU process is proved; this result is extended in [91] to a class of diffusion 
processes admitting steady-state density 

W{x) = hm fix, i |xo ) = ^ exp ( C ^dy) , (96) 

where c is a normalization constant. When the boundary S approaches the 
unattainable level r of the diffusion interval, lima-^^ o'^(a^) [^^(2;)]^ ^{T) = 0, 
the following asymptotic result for the FPT pdf g{t) holds: 

.W-^exp{-^}. (97) 

Numerical studies on the OU and on the Feller processes show that this 
behavior is attained with a negligible error for quite small values of the bound- 
ary S (i.e. for 5' = 3 if /i = 0, 6* = 1 and cr^ = 1 for the OU process). In 
this asymptotic case the mean FPT, E{T), looses the dependency upon the 
initial value xq. Asymptotic results hold for the same processes in the case 
of boundaries cither asymptotically constant or asymptotically periodic |47] . 
Periodic boundaries for the OU process are considered also in 11141 (see [107] 
for a review on time dependent boundaries). 

Let us now switch to the asymptotic behavior with respect to time. For 
small times, the FPT can be approximated with the IG distribution. Indeed 
near the origin any diffusion process can be approximated by a Wiener pro- 
cess. In [47] the asymptotic behaviour, for large t, of the FPT pdf 's through 
some time-varying boundaries is considered. This paper deals with a class of 
one dimensional diffusion processes with steady-state density. The considered 
boundaries include periodic boundaries. Sufficient conditions for an asymp- 
totic exponential behavior are given for the cases of asymptotically constant 
and asymptotically periodic boundaries. Explicit expressions are worked out 
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for the processes that can be obtained from the OU process by spatial trans- 
formations. 

The FPT pdf as i — > oo, for periodic or asymptotically periodic boundaries 
S{t), under mild conditions (cf. [47]) exhibits damped oscillations with the 
same period T as the boundary: 

5(5(t), t\xo) « ait) exp I" ^ a{T)dT^ . (98) 

Here a{t) is a periodic function of period T: 



a(<) = -2 lim i}[S{t + nT),t + nT\xo] (99) 

n— ^oo 

= ~{v'it) + ,[vit)]-j^^^]w[vit)] 

and V{t) = lim„_^oo S{t + nT). This behavior arises also for times not too far 
from the origin. An exponential asymptotic behavior is also proved for large 
times and constant boundaries in |127) . 

An asymptotic evaluation of the probability that the Wiener process first 
crosses a square root boundary is provided in |126j . Denoting as Tc the FPT 
of the Wiener process trough the boundary cy/1 + t one has 

P{T, > t) ^t^eo qt-P^''\ < p(c) < i. (100) 

Here limc^oo p(c) — 0; limc^op(c) = \ and 1p{c) is a real solution between 
and 1 with respect to A of the equation: 

- sin(i^)r(l + A)(V2c)" n-A _ 
^ nnl 2 

n=l 

Using the inverse of transformation (|94l) . this result can be applied to get 
the asymptotic OU process FPT pdf trough a constant boundary for large 
times. 

In [?] truncations of the series expansion of the FPT pdf solution of (|86p 
are used to achieve approximate evaluations. Use of fixed point theorems is 
made to obtain bounds for the FPT pdf of the OU and the Feller processes. 
Inequalities are proved to find for which times the FPT pdf can be approxi- 
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mated, within a preassigned error, by means of an assigned distribution such 
as the FPT of the Wiener process or the exponential one. 

In [93] the asymptotic behavior of the FPT pdf through time-varying 
boundaries is determined for a class of Gauss-Markov processes. 



4.1.4 Moments of the FPT 

Analytical formulae for the moments of the FPT are available only for time 
homogeneous processes with constant boundary. Three approaches are pos- 
sible: 

1. derivatives of the Laplace transform of the FPT pdf; 

2. solution of second order differential equations; 

3. solution of the recursive formula proposed by Siegert ( |130) ). 
Having the Laplace transform of the FPT pdf, one can compute: 

ET- ^ (-1)" (102) 

where g\{S\xQ) is given by (|80|) . The presence of special functions in the 
Laplace transforms for the OU or for the Feller processes (|55l) leads to 
very complex computations. 

Alternatively, using the Kolmogorov equation and eq. ()80|) . one can show 
that the moments of the FPT verify the recursive system of ordinary differ- 
ential equations: 



dxg dxo 



a\xo) + ^i{xa) -r^^ = -nET'^-'ixo), xo G {I, S) (103) 



with boundary conditions: 

ET°{1) = 0, ET°{S) = 1. (104) 

When the process admits steady state distribution, one can write the so- 
lution of (|102[) through the Siegert formula (cf. [130) ): 



ET" = t,,{S\xo)=n r J^f' rW{y)t„MS\y)dy. (105) 
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Due to the numerical difBculties of these formulae, in [84 approximations 
are proposed for specific processes (cf. Subsections 13.51 and I3.6P together with 
suggestions on the use of each one for specific ranges of the parameters. 



4.1.5 Stochastic ordering 

A further technique for the study of the FPT's is the stochastic comparison 
of the FPT's from different models (cf. [123], [118], [Tig). Let us consider 
the FPT's Ti and T2 of two diffusion processes Xi and X2 over Ii = (/i,ri) 
and I2 — (^2i''2) with drifts fj.i{x), x = 1,2 and diffusion coefficients (Ti{x), 
i = 1,2 respectively. Let the two processes Yi and Y2 be obtained from Xi 
and X2 through the transformation 

dz 

y,=gi{x)= ——-,1 = 1,2. 
Moreover, let Yi and Y2 verify the inequalities: 



> ^^YM Vy e [0,32(5)]; < Vy e [0,52(5)] (106) 

ay 

and (jfix) > cr|(a;). 

For xq e (max(/i, Z2), 5), S € (max(Zi. /2), min(ri, r2)), a;o < S, it holds: 



TxiiS\xo) <as Tx^iSlxo). 



(107) 



In (|107p , " as" means " almost surely" . Note that Yi and I2 are characterized 
by unit diffusion coefficien and drift 



fJ-Y, (y) 



1 , 1 daf , , , 



ai{x') 4 dx 



(108) 



4.1.6 Jump diffusion processes 

The following integral equation for the FET pdf 5 ( 5, 1 1 y, r) of the process X 
in ([68]). defined over / = [Z, 5] and originated in the state y at time r, holds 

(cf. my- 
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g{S,t\y,T) 




(109) 



X {A^r {z~a,u\y,T) + X'T {z + a,u \y,T) I^i^s-a) {z)} 



X 5(5,t|z,M) + A'=e-^(*-^) / dzriz,t\y,T). 



S-a 



Here A = A'^ + A% Ia{-) is the indicator function of the set A, the jump 
amphtudes are 5^ = —(5' = a. Furthermore g {S,t\y,T) and /°(x,i|y, s) are 
the FPT and the transition pdf in the presence of the boundary S of the 
underlying diffusion process. The fohowing approximate solution: 



holds (cf. ^9\) for a Wiener process with drift fj, and diffusion coefficient cr 

when A*^ > A*, A*^ ^ 1. Here gt{z, u) = g{^, t\z, u). This approximation can be 
interpreted in terms of sample path behavior for the process X . For jumps of 
low frequency, but relevant amplitude with respect to S, most of the sample 
paths cross the boundary either for diffusion without jumps or for an upward 
jump when Xt (z [S — a, S) or for diffusion after at most a single (upward 
of downward) jump. The possible occurrence of a higher number of jumps 
is disregarded. Hence this approximation explains the first two peaks of the 
observed multimodal behavior exhibited by the FET pdf (cf. Fig. |4]). 

Some results on the moments of two simplified jump diffusion neuronal 
models are discussed in [50] and [51]. In the "large jumps" model the am- 
plitude of exponentially time distributed jumps is large enough to determine 
a crossing of the threshold at each jump. Assuming that the crossing is a 
certain event, a recursive relation holds for the FET moments of order n > 1 
of this model: 
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E[T'^-^] + (-1) 



(111) 



Here gx+g{S \xq) is the Laplace transform, of parameter X + 9, oi the pure 
diffusion FPT, A the frequency of jumps and the superscript (m) denotes the 
m-th derivative with respect to 9. 

In the 'reset model' exponentially time distributed jumps force the mem- 
brane potential to return instantaneously to its resting level Vr = xq. Then 
the dynamics restarts anew till the crossing of the boundary or a new reset- 
ting. This model includes both upward and downward jumps of frequencies 
Ai and A2, whose time epochs are described by means of two independent 
Poisson processes. One has: 



E\T] 



I- g\{S\xo) ^ 
Xg\{S\xQ) 
2E[T] 



(112) 



dgx+eiS\xo] 



gxiS\xo) [gx{S\xo)]' 
with the notations as in (|llip but A = Ai + A2 



d9 



9=0 



4-2 Numerical methods 
4.2.1 The direct FPT problem 

The direct FPT problem requires to determine the FPT pdf for a given 
process assuming the transition pdf and the boundary shape to be known. A 
large literature exists on numerical methods to solve the integral equations 
for the FPT pdf even for time not homogeneous diffusion processes (cf. [7], 
[H] . [55] . [31], [52], [inS], [ini])- The one proposed in [B] seems to be the 
fastest and most efficient. It consists in discretizing when the function 
k {t) is chosen to get a regular kernel for the second kind Volterra equation 
(cf. ([88]) or ([89)1 for the OU and the Feller processes, respectively). Setting 
t ~ to + kh, k = 1,2, h > 0, the discretized solution of eq. ([86]) is 
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g{S {to + h),to + h \xo, to ) - -2^ {S {to + h),to + h \xo, to ) ; (113) 

g{S {to + kh) ,to + kh \xo,to) -2?/; {S {to + kh) ,to + kh |xo,to ) 

fe-i 

+2h J^aiS {to +jh),to+ jh \xo,to) 
i=i 

X V (5* {to + kh) ,to + kh\S {to + jh) , to + jh) 
k = 2,3,... 

A suitable discretization step is necessary to make the numerical integra- 
tion reliable. The numerical algorithm (I114p uses previous integration steps 
to determine the successive values, hence it is important to get good eval- 
uations on the first intervals. A heuristic rule is to execute at least twenty 
integration steps before the maximum of the FPT pdf occurs. In Fig. [71 panel 
A, we show the FPT pdf of a standard OU process with different integration 
steps. The error in the FPT pdf due to a wrong choice for h is enlightened 
in the evaluation of the FPT distribution (Panel B). 

In |113| a strategy is proposed to solve numerically eq. ([55)1 with an appro- 
priate choice of the integration step. To this aim a time-dependent function, 
the FPT Location function, is introduced. 




O 5 10 15 20 O 5 10 15 20 

t (ms) t (ms) 



Fig. 7 FPT pdf (panel A) and FPT probability distribution (panel B) of an OU process 
obtained numerically with fj, = mVms~^ , = 1 mV^ms~^ , = 1 ms~^, 5 = 1 mV; 
integration steps h = 0.045 {continuous line), h = 0.6 {dashed line) 
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4.2.2 The inverse FPT problem 

The inverse FPT problem requires to determine the expression for the bound- 
ary S{t) when the FPT pdf for a diffusion process is known, either in exact 
form or through a sample of FPT's. Two numerical algorithms are proposed 
to solve this problem for the Wiener process in |139| and extended to the OU 
process in |121) . The first algorithm proposes a Monte Carlo procedure to 
approximate the unknown boundary for the Wiener process stepwise. This 
algorithm is reliable and easily implemented but it is computationally ex- 
pensive. The second approach, purely numerical, is computationally more 
attractive and extensions to processes different from the Wiener and the OU 
ones should hold. We present it here in the case of the Wiener process. By 
integrating Fortet's equation ([79]) in x from S{t) to infinity (cf. |97j) one 
obtains the integral equation 

where fFix) = 1 — <f (a;), <P{x) = dz, (^(2) is the standard normal pdf 

and g{t) — g{S{t),t \xo) is the FPT pdf for the Wiener process through the 
threshold S{t). Equation ()114p is a Volterra integral equation of the first kind 
in g(s) but it is a non-linear Volterra integral equation of the second kind 
in S{t). Its kernel is nonsingular since it is bounded. Moreover the functions 
involved in the equation are bounded and 'F is invertible. Hence one can 
obtain numerically the approximate value S*{ti) of S{ti) at the knots ti = ih 
for I = 1, . . . , n; here h = t/n {t > fixed). The integral on the l.h.s. of ()114|) 
is approximated by the Euler method: 



getting a non- linear system of n equations in n unknowns S* (ti), . . . , S* (i„). 
Note that the i-th equation, i >2, makes use of the values S*{ti),. . . ,S*{ti^i) 
in the preceding steps. Hence (|115|) can be solved iteratively to get approxi- 
mate values for the unknown boundary S at the knots. The convergence of 
this algorithm and an estimate of its error are considered in [139] . 
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Extensions to the case in which the FPT pdf is known only through sam- 
ples of ISI's make use of the kernel method to approximate the FPT pdf (cf. 
|122| ). Applications to neuronal modeling are proposed in |121) and |125j for 
the OU process. In particular in |125) this algorithm is employed to propose 
a classification method of groups of neurons when simultaneously recorded 
spike trains are available. 



4-3 Simulation methods 

In the study of neuronal models, when the process is not time homogeneous 
or the boundary is time dependent, the only available general technique is 
simulation. Despite its large use, this methods hides problems that may make 
the results for FPT's unreliable. The standard approach to the simulation of 
FPT's makes use of discretization algorithms for the SDE describing the 
membrane potential dynamics. Various reliable discretization schemes exist 
(cf. |74) and references quoted therein) , depending on which degree of strong 
or of weak convergence is required. The easiest of these schemes is the so 
called Euler-Maruyama one (cf. SUSANNE). 

The major cause of error in the simulation of FPT's is related to the chance 
to leave the crossing of the boundary undetected due to the discretization 
of the sample paths. Indeed a crossing happened inside the discretization 
interval results hidden in the discretized sample (cf. Fig. [8]). This implies 
an important overestimation of the FPT, that does not disappear when the 
discretization step decreases (cf. [53]). The number of trials where the error 
may occur increases with the decrement of h, balancing the corresponding 
decrease in the probability of error. 

Different solutions have been proposed to make the simulation of FPT's 
reliable. They suggest to evaluate the crossing probability during each inte- 
gration step through the bridge process joining the last two values generated 
for the diffusion process. A bridge process ~ t e [io,^i]| is 

associated to a given diffusion X by constraining X to take fixed values at 
the time instants ta and ti > to- The process xl*"'*!! is still a diffusion, since 
its sample paths are a subset of the set of sample paths of X. 
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12 




t (ms) 

Fig. 8 A sample path of a diffusion process and its discretization. Exemplification of a 
missed crossing of the boundary inside a time interval of the simulation. Circles: simulated 
values of the sample path; dots: sample path 

The crossing probability of the original diffusion during a simulation step, 
of amplitude h, coincides with that of its associated bridge on the same time 
interval (cf. for example |112j ). Then one can evaluate, on each interval, the 
probability of hidden crossing for this process. For the Wiener process, the 
probability that the bridge V^I'^^'^+'^l, originated in the state y at time r and 
constrained to assume the value z at time t + h, crosses the boundary S > y 
during [r, t + h] is 

i 2(S^-Sy-Sz + zy)] , , 

P*=Pw'(5,/i,y,z)=expj ^ ^i. (116) 

One can compare this value with a generated random number \J uniform 
on (0,1) and, \i U < P* , conclude that a crossing has happened in that 
interval. The crossing probability of a Wiener process is used to approximate 
the crossing of the bridge associated to the process of interest in To 
introduce a more precise estimation of the crossing probability for the bridge 
we first recall the relationship between the transition pdf's f {x,t\y,T) and 
ix,t\y,T) of the process X and of its bridge (cf. CT): 
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/I*-*^! (x, t\y,T) = f{x,t \y,r) { '"^'^^ to<r<t<h. (117) 

Denoting with r[*"'*il the FPT of the bridge through the boundary 

S and with gl*":*!] (S,t\xo,tQ) its pdf it holds (of. [55]): 

g[*"-*^l (5, t |xo, ^0 ) ^ g (5, ^ |xo, to ) //^: h<T<t<h. (118) 

/ (z,ii \xQ,ta) 

Integral equations analoguous to ([5T|) and (155)) hold also for the FPT 
pdf of the bridge process. In [53] an approximate value of g[*o.*i] (^S,t\xQ) 
is obtained, by disregarding the integral on the l.h.s of such equations. The 
approximation using eq. (|86p 

5^*-*^! (5, t |a;o, <o ) = -2^^ (5, t |xo, to ) ! (119) 

/ (2:,tl Fo,to) 

produces very good estimates in the case of the OU and of the Feller 
underlying diffusion processes. In [53j the integral of this approximation over 
the discretization interval is used to estimate the probability of a hidden 
crossing inside each interval. In [55] a Monte Carlo algorithm is proposed to 
estimate the crossing probability of the bridge process. A numerical scheme 
is applied to the SDE for the bridge process to generate N samples. If L 
samples cross the boundary, the ratio L/N is used to estimate the crossing 
probability. The SDE for the bridge has drift and diffusion coefhcients: 

M[*-*^](x,t) = ^,{x) + ^M^|-/(z,ii|a;,t); frl*-*^! (:r) = a{x) (120) 

/(z,ti|x,t) ax 

respectively (cf. [51]). Here ^(x) and fT^(a;) are the drift and the diffusion 
coefficient of the original diffusion X. 

For a standardized OU process the coefficients (|120p are (cf. [55]): 

^}l°^^\x,t) = -x+^-^^^^^■ al^^%) = l, (121) 
and for the Feller process are (cf. [55]): 
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pP(«i-t) 



p-2- 



eP(*i-0 _ 1 



2pVa;zeP<'i-'> 
r-(eP(*i-*)-l) 



gp(ti-t) _ I 



r(eP('i-*)-l) 



= V2^. (122) 

Here we employed the notation (|511) and /(yy) denotes the nLodificd Bessel 
function of parameter rj. 

A nested algorithm is proposed in |55) for the numerical solution of the 
SDE for the bridge. This choice avoids to evaluate the drift in t — ti where 
it becomes singular. 

An alternative method to evaluate the hidden crossing probabilities, based 
on large deviation arguments, is proposed in [8 . This method is less precise 
than the previousluy mentioned ones but it does not request the knowledge 
of the transition pdf of the the process X. 

These algorithms can be applied also to jump diffusion processes but when- 
ever a jump falls in between the two nodes i„ and t„+i of the partition, the 
right end of the time interval [tmtyi+i] should be substituted with the epoch 
tn, tn < tn < ^ri+i, of the jump evcnt. To account for the possible hidden 
crossings inside the discretization intervals, the correction algorithm proposed 
in [53] should be employed. 

A novel numerical method for the simulation of FPT has been recently pro- 
posed in [133j . The algorithm makes use of the representation of the stochastic 
process through an expansion using the Haar functions. It takes advantage 
of the dichotomic nature of this development to refine the description of the 
process in intervals where posssible hidden crossings may arise. 

In a recent paper [60j it is remarked that the membrane potential, until 
the spike, evolves in the presence of the boundary. The SDE for the process 
constrained by the boundary, i.e. for the process that has not yet attained the 
boundary till a fixed time, is determined. The SDE for its bridge, conditioned 
to cross the boundary for the first time at its right end, is also determined. Use 
of the simulation techniques allows to simulate these stochastic differential 
equations. 
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A few papers exist on the parameter estimation problem. The hterature on 
this subject is rather recent, disregarding two oldest papers. The first one 
([77]) considers a sample of membrane potential values observed at discrete 
times while [55] uses the moment method on a sample of ISI's. We focus here 
mainly on the available statistical methods for the OU and the Feller models. 
Their parameters can be divided into two groups: the intrinsic parameters, 
S, xq, Vj and 9 for the OU process and S,Xo,t for the Feller process, and the 
input parameters, /i and cr^ for the OU process and and ct| for the Feller 
one. The intrinsic parameters are often disregarded in estimation problems 
assuming their direct measure. In j63j the estimation of the refractory period 
is also discussed. 

We distinguish in the sequel two types of methods, depending upon the 
available data: 

1. Intracellular membrane recordings; 

2. ISI time series. 



5.1 Samples from membrane potential measures 

In |85| a regression method and a maximum likelihood technique are applied 
to estimate /3 = ^, /i and a from intracellular data, supposed to follow an 
OU process. We report here the maximum likelihood estimators for the case 
of OU and Feller processes. One assumes that during an ISI the membrane 
depolarization Xi = Xi, i = 0, 1, N, is sampled at the N + 1 points U = ih. 
The maximum likelihood estimates are: 

^ 22j=o xj+x N 

n-^^^+Px, (124) 

and 
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T 



1 - \^ 

a = (xj+i - Xj + Xjhfi - h^j (125) 



3=0 



where x = Sj=o ^J' T = Nh. The hkehhood estimators 'jlp and a for 
the parameters ct — ^ and j^ip of the FcUer process (of. 12|) coincide with 
(|123p and (|124p . while the estimator for tT|, putting ft/i = ( 1 — e^^ ) , is: 



^2 = \ ^- (126) 

Formulae (jl23p - (|125p are obtained disregarding the existence of the firing 
boundary. This approximation determines a bias on the estimated values (cf. 
[S] , [in] ) • The bias of the estimator of ji is of the same order of magnitude as 
the standard deviation of the estimate. 

Unbiased estimators for /x are not yet available for the OU and the Feller 
models while the bias for the RRW and for the Wiener models are computed 
in a closed form in [5]. A comparative study on the estimators for the Feller 
process is performed in jlO| . In |11| maximum likelihood estimators are de- 
rived from discrete observations of a Markov process up to the first-hitting 
time of a threshold, both in discrete and in continuous time. The models 
considered are the RRW, an autoregressive model of order one (AR(1)), and 
the Wiener, OU and Feller diffusions. For the last two ones approximations 
are introduced to evaluate the conditional transition pdf and the FPT dis- 
tribution. These approximations hold when the sampling intervals are small. 
Their use allow to evaluate the likelihood function. 

In [55] an algorithm is proposed to compute likelihoods, based on the 
numerical solution of the integral equation (j86p . Furthermore an estimator 
based on the large deviation principle is suggested to deal with the case of 
very small likelihoods in the tails of the distribution. 

In |98| a maximum likelihood estimation method is employed for a par- 
ticular LIF model with an additional variance parameter modeling possible 
slow fluctuations in the parameter ^. 

In [67] a sample of discrete observations X^/i, io < i < ii, io '■— \^ '■— 
[3! process X of eq. ([5^ . over the time interval [to,^i], is considered. 

The following nonparametric kernel estimators are proposed: 
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^ i=io V h ) \ AM J 



(72 (a) = ' X ^ (127) 

with /i > a suitable bandwith. Furthermore K (y) may be chosen as 
a rectangular or triangular kernel, for a suitable integer M. Examples of 
possible kernels are K{y) = (y) and K{y) = (1 - |y|) /{-i.+i} (y), 

with /{,}(j/) indicator function of the set {•} . 

Extensions to other processes and for the selection of the model (OU, 
Feller or different ones) are used to exemplify the method. The examples 
considered correspond to rarely spiking neurons, a fact that minimizes the 
problem underlined in [9], but prevents its use in other instances. 



5.2 Samples of ISI's 

The case of ISI data has been recently considered in [90_. In this paper 
an algorithm is proposed for computing maximum likelihood estimates with 
their confidence regions for /i and a^. The algorithm numerically inverts the 
Laplace transform for the OU model. The method works also to estimate the 
parameter 6 but it requests larger samples. 

Maximum likelihood estimates for the OU model are obtained in [138j 
using numerical evaluations of (p3|) . In [32] a variant of the moment method 
is proposed to estimate the input parameters of the OU process. In this paper 
an optimal stopping theorem is applied to determine the first two exponential 
moments of the FPT: 



E 



The moment method is then applied to obtain the estimators /i„ and tr^ 
from a sample of ISI's {Ti, T„}: 
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where 

-J n 1 " 

= i V Z2.„ = i V e^^^ir (130) 

i=\ i=l 

This method can be apphed only in the suprathreshold region since the 
following conditions must be fulfilled: E (e-^/^) < oo, E (e^^^^) < oo. The 
first condition is verified if ^9 > S and the second holds if fi6 > S and 

In [33] analoguous results are proved for the Feller model, for which 



V ^ MFT - S" V ) (^^pr - sf - air (/ij.T/2 - 5) ' 

(131) 

These expectations converge when ^pT > S and ^y^l + — 1^ < 
(yLiii^T — 5). The estimators are: 



Mi^,™ = — ^ r (1^2) 



r [2 (Zi,„ - 1) {SZ2,n - yo) - (5*^1, n - yo) (^2,n - 1)] 

Consistency and asymptotic normality of these estimators have been 
proved in [5S] . The sample sizes required to get the asymptotic conditions are 
not huge (some hundreds), hence this property can be performed in neuronal 
experiments. 

In |34| an alternative method to estimate /l( and a^, based on the analogu- 
ous of (|114p for the OU process, 



/ 



/i6i(l-e-^)-5' 



sll ( 

2 V 



1 — e 6 



aVO V 1 + e- 



du, (133) 



is proposed. Numerical results suggest that this approach is preferable to 
the previous ones. 
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The case of observations of the trajectory on very short time intervals is 
considered in [SS] . They propose a method to estimate the parameters of an 
OU process in this particular instance. 

A recent review ([86]) has appeared summarizing the state of the art of 
the estimation problem for diffusion processes in neuromodeling instances. A 
comparison of the different estimation methods for the OU process is per- 
formed in [55] . 

Parallel spike trains are deeply discussed in |61) . This book presents the 
methods of correlation analysis together with a review on different approaches 
for the analysis of single spike trains. The different chapters discuss many 
important problems related with the statistical analysis of spike trains. 

Finally the inverse FPT method is applied in [125] to classify simultane- 
ously recorded spike trains. The value of the parameters for the OU pro- 
cess are assumed constant for all the recorded spike trains. The boundary 
is determined by the inverse FPT method and a comparison of the different 
boundaries is employed to classify the data. 

The case of not stationary data is not contemplated by the estimation 
procedure. The problem of models whose noise term has a specific temporal 
structure has not been solved up to now. In [122j the inverse FPT method 
is used on samples of FPT's from an OU process to recover the boundary 
shape and to test nonstationary behaviors. The proposed method makes use 
of a moving window approach. The inverse FPT is applied on samples from 
each window. The comparison of the determined boundary shapes allows to 
detect changes in the observed dynamics. 
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